


#######################################################################################

MeanOutFun <- function(IN_FILE, OUT_FILE){

result.array <- readRDS(paste("/Users/richardtraunmuller/Dropbox/rep_hegsamb/", IN_FILE, ".rds", sep=""))
mean.array <- matrix(NA, 6, 8)
rownames(mean.array) <- rownames(result.array)
colnames(mean.array) <- c("NumRob", "TP", "FP", "TN", "FN", "CSign", "RMSE", "Corrs")
for(i in 1:6){
  for(j in 1:8){
    
    mean.array[i, j] <- mean(result.array[i, j, ], na.rm=T) 
    
  }
}

write.table(round(mean.array, 2), file=paste("/Users/richardtraunmuller/Dropbox/PSRM_Replication/", OUT_FILE,  ".csv", sep=""), sep=",")
print(mean.array)
}

###############################################################################
# One Variable
MeanOutFun("one_result_mat_1", "one_means_25")
MeanOutFun("one_result_mat_2", "one_means_50")
MeanOutFun("one_result_mat_3", "one_means_75")
MeanOutFun("one_result_mat_4", "one_means_10")
MeanOutFun("one_result_mat_5", "one_means_125")

# Five Variables
MeanOutFun("five_result_mat_1", "five_means_25")
MeanOutFun("five_result_mat_2", "five_means_50")
MeanOutFun("five_result_mat_3", "five_means_75")
MeanOutFun("five_result_mat_4", "five_means_10")
MeanOutFun("five_result_mat_5", "five_means_125")

# Nine Variables
MeanOutFun("nine_result_mat_1", "nine_means_25")
MeanOutFun("nine_result_mat_2", "nine_means_50")
MeanOutFun("nine_result_mat_3", "nine_means_75")
MeanOutFun("nine_result_mat_4", "nine_means_10")
MeanOutFun("nine_result_mat_5", "nine_means_125")

###############################################################################
# One Runif Variable
MeanOutFun("one_runif_result_mat_1", "one_runif_means_25")
MeanOutFun("one_runif_result_mat_2", "one_runif_means_50")
MeanOutFun("one_runif_result_mat_3", "one_runif_means_75")
MeanOutFun("one_runif_result_mat_4", "one_runif_means_10")
MeanOutFun("one_runif_result_mat_5", "one_runif_means_125")

# Five Runif Variable
MeanOutFun("five_runif_result_mat_1", "five_runif_means_25")
MeanOutFun("five_runif_result_mat_2", "five_runif_means_50")
MeanOutFun("five_runif_result_mat_3", "five_runif_means_75")
MeanOutFun("five_runif_result_mat_4", "five_runif_means_10")
MeanOutFun("five_runif_result_mat_5", "five_runif_means_125")

# Nine Runif Variable
MeanOutFun("nine_runif_result_mat_1", "nine_runif_means_25")
MeanOutFun("nine_runif_result_mat_2", "nine_runif_means_50")
MeanOutFun("nine_runif_result_mat_3", "nine_runif_means_75")
MeanOutFun("nine_runif_result_mat_4", "nine_runif_means_10")
MeanOutFun("nine_runif_result_mat_5", "nine_runif_means_125")

##############################################################################
# Corr 5/10 norm 1
MeanOutFun("five_norm_1_result_mat_1", "five_norm_1_means_1")
MeanOutFun("five_norm_1_result_mat_2", "five_norm_1_means_2")
MeanOutFun("five_norm_1_result_mat_3", "five_norm_1_means_3")
MeanOutFun("five_norm_1_result_mat_4", "five_norm_1_means_4")
MeanOutFun("five_norm_1_result_mat_5", "five_norm_1_means_5")

# Corr 5/10 norm 3
MeanOutFun("five_norm_3_result_mat_1", "five_norm_3_means_1")
MeanOutFun("five_norm_3_result_mat_2", "five_norm_3_means_2")
MeanOutFun("five_norm_3_result_mat_3", "five_norm_3_means_3")
MeanOutFun("five_norm_3_result_mat_4", "five_norm_3_means_4")
MeanOutFun("five_norm_3_result_mat_5", "five_norm_3_means_5")

# Corr 5/10 norm 5
MeanOutFun("five_norm_5_result_mat_1", "five_norm_5_means_1")
MeanOutFun("five_norm_5_result_mat_2", "five_norm_5_means_2")
MeanOutFun("five_norm_5_result_mat_3", "five_norm_5_means_3")
MeanOutFun("five_norm_5_result_mat_4", "five_norm_5_means_4")
MeanOutFun("five_norm_5_result_mat_5", "five_norm_5_means_5")
###############################################################################


# The probability of false negatives is 1 minus 
# the average number of true positives divided by 
# the number of determinants p.

# The probability of false positives is one minus the average number 
# of true negatives divided by the number of confounders k-p. 

################################################################################
# Table 1: Probability of False Positive and False Negatives When Variables are Uncorrelated
p <- 1

n1 <- round(1-MeanOutFun("one_result_mat_1", "one_means_25")[c(1:2, 6, 3:5), c(2)]/p, 3)
n2 <- round(1-MeanOutFun("one_result_mat_2", "one_means_50")[c(1:2, 6, 3:5), c(2)]/p, 3)
n3 <- round(1-MeanOutFun("one_result_mat_3", "one_means_75")[c(1:2, 6, 3:5), c(2)]/p, 3)
n4 <- round(1-MeanOutFun("one_result_mat_4", "one_means_10")[c(1:2, 6, 3:5), c(2)]/p, 3)
n5 <- round(1-MeanOutFun("one_result_mat_5", "one_means_125")[c(1:2, 6, 3:5), c(2)]/p, 3)

p1 <- round(1-MeanOutFun("one_result_mat_1", "one_means_25")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)
p2 <- round(1-MeanOutFun("one_result_mat_2", "one_means_50")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)
p3 <- round(1-MeanOutFun("one_result_mat_3", "one_means_75")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)
p4 <- round(1-MeanOutFun("one_result_mat_4", "one_means_10")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)
p5 <- round(1-MeanOutFun("one_result_mat_5", "one_means_125")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)

n <- t(rbind(n1, n2, n3, n4, n5))
p <- t(rbind(p1, p2, p3, p4, p5))

tab.1.a <- rbind(n[1,], p[1,], n[2,], p[2,], n[3,], p[3,], n[4,], p[4,], n[5,], p[5,], n[6,], p[6,])

#
p <- 5

n1 <- round(1-MeanOutFun("five_result_mat_1", "five_means_25")[c(1:2, 6, 3:5), c(2)]/p, 3)
n2 <- round(1-MeanOutFun("five_result_mat_2", "five_means_50")[c(1:2, 6, 3:5), c(2)]/p, 3)
n3 <- round(1-MeanOutFun("five_result_mat_3", "five_means_75")[c(1:2, 6, 3:5), c(2)]/p, 3)
n4 <- round(1-MeanOutFun("five_result_mat_4", "five_means_10")[c(1:2, 6, 3:5), c(2)]/p, 3)
n5 <- round(1-MeanOutFun("five_result_mat_5", "five_means_125")[c(1:2, 6, 3:5), c(2)]/p, 3)

p1 <- round(1-MeanOutFun("five_result_mat_1", "five_means_25")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)
p2 <- round(1-MeanOutFun("five_result_mat_2", "five_means_50")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)
p3 <- round(1-MeanOutFun("five_result_mat_3", "five_means_75")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)
p4 <- round(1-MeanOutFun("five_result_mat_4", "five_means_10")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)
p5 <- round(1-MeanOutFun("five_result_mat_5", "five_means_125")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)

n <- t(rbind(n1, n2, n3, n4, n5))
p <- t(rbind(p1, p2, p3, p4, p5))

tab.1.b <- rbind(n[1,], p[1,], n[2,], p[2,], n[3,], p[3,], n[4,], p[4,], n[5,], p[5,], n[6,], p[6,])

#
p <- 9

n1 <- round(1-MeanOutFun("nine_result_mat_1", "nine_means_25")[c(1:2, 6, 3:5), c(2)]/p, 3)
n2 <- round(1-MeanOutFun("nine_result_mat_2", "nine_means_50")[c(1:2, 6, 3:5), c(2)]/p, 3)
n3 <- round(1-MeanOutFun("nine_result_mat_3", "nine_means_75")[c(1:2, 6, 3:5), c(2)]/p, 3)
n4 <- round(1-MeanOutFun("nine_result_mat_4", "nine_means_10")[c(1:2, 6, 3:5), c(2)]/p, 3)
n5 <- round(1-MeanOutFun("nine_result_mat_5", "nine_means_125")[c(1:2, 6, 3:5), c(2)]/p, 3)

p1 <- round(1-MeanOutFun("nine_result_mat_1", "nine_means_25")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)
p2 <- round(1-MeanOutFun("nine_result_mat_2", "nine_means_50")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)
p3 <- round(1-MeanOutFun("nine_result_mat_3", "nine_means_75")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)
p4 <- round(1-MeanOutFun("nine_result_mat_4", "nine_means_10")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)
p5 <- round(1-MeanOutFun("nine_result_mat_5", "nine_means_125")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)

n <- t(rbind(n1, n2, n3, n4, n5))
p <- t(rbind(p1, p2, p3, p4, p5))

tab.1.c <- rbind(n[1,], p[1,], n[2,], p[2,], n[3,], p[3,], n[4,], p[4,], n[5,], p[5,], n[6,], p[6,])

table.1 <- rbind(tab.1.a, tab.1.b, tab.1.c)
rownames(table.1) <- rep(rep(c("EBA_Leamer", "EBA_sign", "EBA_90", "MA_unig", "MA_Sala", "MA_BIC"), each=2), 3)
colnames(table.1) <- c(".25", ".50", ".75", "1.00", "1.25")

write.table(table.1, file="/Users/richardtraunmuller/Dropbox/PSRM_Replication/pluemper_traunmueller_table1.csv", sep=",")

##################################################################################################################################
# Table 2: Probability of False Positive and False Negatives When Variables are Highly Correlated
p <- 1

n1 <- round(1-MeanOutFun("one_runif_result_mat_1", "one_runif_means_25")[c(1:2, 6, 3:5), c(2)]/p, 3)
n2 <- round(1-MeanOutFun("one_runif_result_mat_2", "one_runif_means_50")[c(1:2, 6, 3:5), c(2)]/p, 3)
n3 <- round(1-MeanOutFun("one_runif_result_mat_3", "one_runif_means_75")[c(1:2, 6, 3:5), c(2)]/p, 3)
n4 <- round(1-MeanOutFun("one_runif_result_mat_4", "one_runif_means_10")[c(1:2, 6, 3:5), c(2)]/p, 3)
n5 <- round(1-MeanOutFun("one_runif_result_mat_5", "one_runif_means_125")[c(1:2, 6, 3:5), c(2)]/p, 3)

p1 <- round(1-MeanOutFun("one_runif_result_mat_1", "one_runif_means_25")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)
p2 <- round(1-MeanOutFun("one_runif_result_mat_2", "one_runif_means_50")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)
p3 <- round(1-MeanOutFun("one_runif_result_mat_3", "one_runif_means_75")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)
p4 <- round(1-MeanOutFun("one_runif_result_mat_4", "one_runif_means_10")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)
p5 <- round(1-MeanOutFun("one_runif_result_mat_5", "one_runif_means_125")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)

n <- t(rbind(n1, n2, n3, n4, n5))
p <- t(rbind(p1, p2, p3, p4, p5))

tab.2.a <- rbind(n[1,], p[1,], n[2,], p[2,], n[3,], p[3,], n[4,], p[4,], n[5,], p[5,], n[6,], p[6,])

#
p <- 5

n1 <- round(1-MeanOutFun("five_runif_result_mat_1", "five_runif_means_25")[c(1:2, 6, 3:5), c(2)]/p, 3)
n2 <- round(1-MeanOutFun("five_runif_result_mat_2", "five_runif_means_50")[c(1:2, 6, 3:5), c(2)]/p, 3)
n3 <- round(1-MeanOutFun("five_runif_result_mat_3", "five_runif_means_75")[c(1:2, 6, 3:5), c(2)]/p, 3)
n4 <- round(1-MeanOutFun("five_runif_result_mat_4", "five_runif_means_10")[c(1:2, 6, 3:5), c(2)]/p, 3)
n5 <- round(1-MeanOutFun("five_runif_result_mat_5", "five_runif_means_125")[c(1:2, 6, 3:5), c(2)]/p, 3)

p1 <- round(1-MeanOutFun("five_runif_result_mat_1", "five_runif_means_25")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)
p2 <- round(1-MeanOutFun("five_runif_result_mat_2", "five_runif_means_50")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)
p3 <- round(1-MeanOutFun("five_runif_result_mat_3", "five_runif_means_75")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)
p4 <- round(1-MeanOutFun("five_runif_result_mat_4", "five_runif_means_10")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)
p5 <- round(1-MeanOutFun("five_runif_result_mat_5", "five_runif_means_125")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)

n <- t(rbind(n1, n2, n3, n4, n5))
p <- t(rbind(p1, p2, p3, p4, p5))

tab.2.b <- rbind(n[1,], p[1,], n[2,], p[2,], n[3,], p[3,], n[4,], p[4,], n[5,], p[5,], n[6,], p[6,])

#
p <- 9

n1 <- round(1-MeanOutFun("nine_runif_result_mat_1", "nine_runif_means_25")[c(1:2, 6, 3:5), c(2)]/p, 3)
n2 <- round(1-MeanOutFun("nine_runif_result_mat_2", "nine_runif_means_50")[c(1:2, 6, 3:5), c(2)]/p, 3)
n3 <- round(1-MeanOutFun("nine_runif_result_mat_3", "nine_runif_means_75")[c(1:2, 6, 3:5), c(2)]/p, 3)
n4 <- round(1-MeanOutFun("nine_runif_result_mat_4", "nine_runif_means_10")[c(1:2, 6, 3:5), c(2)]/p, 3)
n5 <- round(1-MeanOutFun("nine_runif_result_mat_5", "nine_runif_means_125")[c(1:2, 6, 3:5), c(2)]/p, 3)

p1 <- round(1-MeanOutFun("nine_runif_result_mat_1", "nine_runif_means_25")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)
p2 <- round(1-MeanOutFun("nine_runif_result_mat_2", "nine_runif_means_50")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)
p3 <- round(1-MeanOutFun("nine_runif_result_mat_3", "nine_runif_means_75")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)
p4 <- round(1-MeanOutFun("nine_runif_result_mat_4", "nine_runif_means_10")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)
p5 <- round(1-MeanOutFun("nine_runif_result_mat_5", "nine_runif_means_125")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)

n <- t(rbind(n1, n2, n3, n4, n5))
p <- t(rbind(p1, p2, p3, p4, p5))

tab.2.c <- rbind(n[1,], p[1,], n[2,], p[2,], n[3,], p[3,], n[4,], p[4,], n[5,], p[5,], n[6,], p[6,])

table.2 <- rbind(tab.2.a, tab.2.b, tab.2.c)
rownames(table.2) <- rep(rep(c("EBA_Leamer", "EBA_sign", "EBA_90", "MA_unig", "MA_Sala", "MA_BIC"), each=2), 3)
colnames(table.2) <- c(".25", ".50", ".75", "1.00", "1.25")

write.table(table.2, file="/Users/richardtraunmuller/Dropbox/PSRM_Replication/pluemper_traunmueller_table2.csv", sep=",")

##################################################################################################################################
# Table 3: Probability of False Positive and False Negatives for Varying Levels of Correlation 
p <- 5

n1 <- round(1-MeanOutFun("five_norm_1_result_mat_1", "five_norm_1_means_25")[c(1:2, 6, 3:5), c(2)]/p, 3)
n2 <- round(1-MeanOutFun("five_norm_1_result_mat_2", "five_norm_1_means_50")[c(1:2, 6, 3:5), c(2)]/p, 3)
n3 <- round(1-MeanOutFun("five_norm_1_result_mat_3", "five_norm_1_means_75")[c(1:2, 6, 3:5), c(2)]/p, 3)
n4 <- round(1-MeanOutFun("five_norm_1_result_mat_4", "five_norm_1_means_10")[c(1:2, 6, 3:5), c(2)]/p, 3)
n5 <- round(1-MeanOutFun("five_norm_1_result_mat_5", "five_norm_1_means_125")[c(1:2, 6, 3:5), c(2)]/p, 3)

p1 <- round(1-MeanOutFun("five_norm_1_result_mat_1", "five_norm_1_means_25")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)
p2 <- round(1-MeanOutFun("five_norm_1_result_mat_2", "five_norm_1_means_50")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)
p3 <- round(1-MeanOutFun("five_norm_1_result_mat_3", "five_norm_1_means_75")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)
p4 <- round(1-MeanOutFun("five_norm_1_result_mat_4", "five_norm_1_means_10")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)
p5 <- round(1-MeanOutFun("five_norm_1_result_mat_5", "five_norm_1_means_125")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)

n <- t(rbind(n1, n2, n3, n4, n5))
p <- t(rbind(p1, p2, p3, p4, p5))

tab.3.a <- rbind(n[1,], p[1,], n[2,], p[2,], n[3,], p[3,], n[4,], p[4,], n[5,], p[5,], n[6,], p[6,])

#
p <- 5

n1 <- round(1-MeanOutFun("five_norm_3_result_mat_1", "five_norm_3_means_25")[c(1:2, 6, 3:5), c(2)]/p, 3)
n2 <- round(1-MeanOutFun("five_norm_3_result_mat_2", "five_norm_3_means_50")[c(1:2, 6, 3:5), c(2)]/p, 3)
n3 <- round(1-MeanOutFun("five_norm_3_result_mat_3", "five_norm_3_means_75")[c(1:2, 6, 3:5), c(2)]/p, 3)
n4 <- round(1-MeanOutFun("five_norm_3_result_mat_4", "five_norm_3_means_10")[c(1:2, 6, 3:5), c(2)]/p, 3)
n5 <- round(1-MeanOutFun("five_norm_3_result_mat_5", "five_norm_3_means_125")[c(1:2, 6, 3:5), c(2)]/p, 3)

p1 <- round(1-MeanOutFun("five_norm_3_result_mat_1", "five_norm_3_means_25")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)
p2 <- round(1-MeanOutFun("five_norm_3_result_mat_2", "five_norm_3_means_50")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)
p3 <- round(1-MeanOutFun("five_norm_3_result_mat_3", "five_norm_3_means_75")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)
p4 <- round(1-MeanOutFun("five_norm_3_result_mat_4", "five_norm_3_means_10")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)
p5 <- round(1-MeanOutFun("five_norm_3_result_mat_5", "five_norm_3_means_125")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)

n <- t(rbind(n1, n2, n3, n4, n5))
p <- t(rbind(p1, p2, p3, p4, p5))

tab.3.b <- rbind(n[1,], p[1,], n[2,], p[2,], n[3,], p[3,], n[4,], p[4,], n[5,], p[5,], n[6,], p[6,])

#
p <- 5

n1 <- round(1-MeanOutFun("five_norm_5_result_mat_1", "five_norm_5_means_25")[c(1:2, 6, 3:5), c(2)]/p, 3)
n2 <- round(1-MeanOutFun("five_norm_5_result_mat_2", "five_norm_5_means_50")[c(1:2, 6, 3:5), c(2)]/p, 3)
n3 <- round(1-MeanOutFun("five_norm_5_result_mat_3", "five_norm_5_means_75")[c(1:2, 6, 3:5), c(2)]/p, 3)
n4 <- round(1-MeanOutFun("five_norm_5_result_mat_4", "five_norm_5_means_10")[c(1:2, 6, 3:5), c(2)]/p, 3)
n5 <- round(1-MeanOutFun("five_norm_5_result_mat_5", "five_norm_5_means_125")[c(1:2, 6, 3:5), c(2)]/p, 3)

p1 <- round(1-MeanOutFun("five_norm_5_result_mat_1", "five_norm_5_means_25")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)
p2 <- round(1-MeanOutFun("five_norm_5_result_mat_2", "five_norm_5_means_50")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)
p3 <- round(1-MeanOutFun("five_norm_5_result_mat_3", "five_norm_5_means_75")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)
p4 <- round(1-MeanOutFun("five_norm_5_result_mat_4", "five_norm_5_means_10")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)
p5 <- round(1-MeanOutFun("five_norm_5_result_mat_5", "five_norm_5_means_125")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)

n <- t(rbind(n1, n2, n3, n4, n5))
p <- t(rbind(p1, p2, p3, p4, p5))

tab.3.c <- rbind(n[1,], p[1,], n[2,], p[2,], n[3,], p[3,], n[4,], p[4,], n[5,], p[5,], n[6,], p[6,])

table.3 <- rbind(tab.3.a, tab.3.b, tab.3.c)
rownames(table.3) <- rep(rep(c("EBA_Leamer", "EBA_sign", "EBA_90", "MA_unig", "MA_Sala", "MA_BIC"), each=2), 3)
colnames(table.3) <- c(".25", ".50", ".75", "1.00", "1.25")

write.table(table.3, file="/Users/richardtraunmuller/Dropbox/PSRM_Replication/pluemper_traunmueller_table3.csv", sep=",")

############################################################################################################################


